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Abstract 

We present a systematic and comprehensive study of finite-size effects in diffusion quantum Monte 
Carlo calculations of metals. Several previously introduced schemes for correcting finite-size errors 
are compared for accuracy and efficiency and practical improvements are introduced. In particular, 
we test a simple but efficient method of finite-size correction based on an accurate combination 
of twist averaging and density functional theory. Our diffusion quantum Monte Carlo results for 
lithium and aluminum, as examples of metallic systems, demonstrate excellent agreement between 
all of the approaches considered. 
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I. INTRODUCTION 


Density functional theory (DFT) has dominated atomic-scale materials modeling for the 
past three decades and will continue to be extremely important. However, there are plenty 
of problems where the accuracy of DFT falls short of requirements. There are clear reasons 
to expect DFT to struggle in strongly-correlated solids, but these are not the only systems 
for which DFT is insufficient. Take, for example, the problem of distinguishing between 
molecular crystal phases and competing low-energy polymorphs. Even in relatively simple 
molecular solids such as crystalline benzene and its polymorphs under pressure, the energy 
differences of interest are less than a few kJ/mol. The most successful calculations based on 
DFT are only reliable to ~10 kJ/mol,^and it has recently been shown that the use of ab initio 
many-electron wavefunction methods, such as quantum Monte Carlo (QMC), is essential to 
tackle this problem successfully.^ Unlike DFT, many-electron wavefunction methods can, in 
principle, be improved systematically until the required convergence is obtained. 

Another important example is the adsorption of molecules on surfaces, where DFT is 
sometimes unable to give predictions of useful accuracy. DFT values of surface formation en¬ 
ergies of simple paradigmatic materials such as silicon and magnesium oxide depend strongly 
on the assumed exchange-correlation functional, and there is usually no way of knowing in 
advance which functional to trust. The well known problem of calculating electronic band 
gaps could also be mentioned. 

The urgent practical need to go beyond DFT in these and other areas is driving current 
efforts to develop more accurate methods. There is abundant evidence that there are large 
classes of problems for which QMC techniques, in particular diffusion Monte Carlo (DMC), 
are considerably more accurate that DFT.®^^ Just recently, using QMC methods, chemically 
accurate ionization potentials have been obtained for the hrst-row transition-metal atoms 
with a mean absolute error of only 0.126 kcal/mol.^By combining high accuracy with good 
efficiency and scalability, QMC methods promise to bring high accuracy to computational 
materials science as a matter of routine.^ 

QMC simulations of extended systems are carried out using hnite simulation cells subject 
to periodic boundary conditions. Practical and computational constraints restrict the size 
of the simulation cell and so introduce hnite-size (FS) errors, which can be large and are 
one of the main problems holding back the application of accurate QMC techniques to 


2 


solids.®^ Quantifying and correcting these errors is an essential part of all QMC simulations 
of extended systems, particularly when high accuracy is required. 

FS errors affect independent-particle approaches such as DFT as well as many-body 
approaches such as QMC. For calculations of perfect crystals, DFT FS errors can be reduced 
simply by improving the accuracy of the Brillouin zone integration, although the errors arise 
when periodic supercells are used to model aperiodic systems are less easily dealt with. 
Independent-particle FS errors also affect many-body calculations of perfect crystalJ^ and 
can be reduced using twist averaging, which is the many-electron equivalent of Brillouin 
zone integration. Even after twist averaging, however, “many-body” FS errors with no 
independent-particle analogues remain. 

Another important contribution to the FS error in many-body simulations of extended 
systems arises from the treatment of the potential energy. The 1/r Coulomb interaction 
is inconsistent with the periodicity of the simulation cell and has to be replaced by the 
Ewald interaction, which is the Green’s function of Poisson’s equation subject to periodic 
boundary conditions. Unlike the Coulomb interaction, the Ewald interaction depends on the 
size and shape of the simulation cell, leading to additional hnite-size errors.^ One approach 
to circumventing this problem is to use a different periodic function, the “model periodic 
Coulomb” (MPC) interaction, in place of the Ewald interaction.^^IIS] Variational quantum 
Monte Carlo (VMC) simulations using the MPC interaction suffer from smaller FS effects 
than simulations using the standard Ewald interaction. It has also been shown that using the 
MPC interaction reduces the FS errors in DMC calculations of ground and excited states.^ 

A drawback of the MPC approach is that it only reduces FS errors arising from the 
use of the Ewald interaction. The charge density and exchange-correlation hole in a hnite 
simulation cell are often very similar to those of an inhnite solid, so the errors in the Coulomb 
energy are indeed primarily due to the errors in the interaction; but the imposition of periodic 
boundary conditions also affects the many-electron wavefunction and thus the electronic 
kinetic energy. The non-interacting part of the FS error in the kinetic energy may be 
eliminated by twist averaging, but the many-body contributions are not negligible. 

Under the assumption that the low-/c behaviour of the structure factor is independent 
of the choice of simulation cell, Chiesa et al^ proposed a method to estimate the many- 
body contributions to the FS errors in both the potential and kinetic energies without 
abandoning the Ewald interaction. Employing this correction, which is based on the random- 
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phase approximation at long wavelength, one can calculate FS corrections within a single 
simulation. 

Another approach to the treatment of many-body hnite-size errors is provided by the 
Kwee-Zhang-Krakauer (KZK) functional,^ which adds a correction computed from the dif¬ 
ference between the DFT energy evaluated using the local density approximation (LDA) for 
an inhnite system and the DFT energy evaluated using a modihed LDA specihcally designed 
to reproduce the total energy of the hnite simulation cell, including FS errors. Both the 
standard LDA and the KZK LDA are parameterized on the basis of DMC simulations of 
cells of uniform electron gas subject to periodic boundary conditions, but the standard LDA 
uses DMC energies that have been extrapolated to inhnite cell size while the KZK LDA 
does not. 

In this paper, we systematically study the problem of eliminating FS errors from QMC 
calculations of real metallic systems, taking lithium and aluminum as examples. We an¬ 
alyze twist-averaged DMC energies obtained using the Ewald interaction and the MFC 
interaction j^^l^^l and hnite-size corrections based on the Chiesa formalismS^^ and the KZK 
functional.^ We also investigate DFT-based corrections designed to improve imperfectly 
twist-averaged results and consider how best to combine the use of twist-averaged boundary 
conditiond^ with the KZK functional.® 


II. COMPUTATIONAL DETAILS 

A. Diffusion quantum Monte Carlo calculations 

The diffusion quantum Monte Carlo method is a stochastic technique for obtaining the 
ground-state energy of a many-electron system. DMC has been described in many previous 
paper^^® and will not be discussed in detail here, but since this work is focused on technical 
aspects of DMC simulations we start with a brief explanation. 

The DMC algorithm solves the imaginary-time Schrodinger equation, 

= I E - (^(R) - t), (1) 

i=l 

where R = (ri, r 2 ,..., r^r^) is a SW-dimensional vector dehning the positions of all W 
electrons in the simulation cell, r is the imaginary time (a real variable despite its name). 
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V (R) is the potential energy including electron-electron interactions, and Et is a constant 
energy offset. (We work in Hartree atomic units, where the numerical values of h, e, rrie, 
and dvr^o are all equal to 1.) The imaginary-time Schrodinger equation resembles a 3Ng- 
dimensional diffusion equation with diffusion constant D = 1/2. The potential energy 


term causes the diffusers to “branch” (multiply or die out) at a position-dependent rate. 


The wavefunction T (R, r) is the number density of diffusers, which are normally known 
as walkers or configurations and are points in the 3W-dimensional configuration space, not 
individual electrons. The DMC algorithm uses this simple physical interpretation to simulate 
the imaginary-time evolution of the wavefunction using a finite population of diffusing and 
branching walkers. 

Solving the imaginary-time Schrodinger equation is useful because it projects out the 
ground state as r —)■ oo. If the initial wavefunction is expanded as a linear combination of 
energy eigenfunctions, T(r = 0) = solution of the imaginary-time Schrodinger 

equation /dr = —{H — is 



( 2 ) 


Thus, as long as Cq ^ 0, the wavefunction T(r) becomes proportional to Tq as r —)■ oo. By 
gradually adjusting Et to maintain the normalization of the solution in the large r limit, 
we can find the ground-state energy Eq. 

An obvious difficulty with this approach is that the wavefunction \k(R, r), which is not 
necessarily positive, is interpreted as a walker density, which must be positive. In fact, a 
naive application of the DMC algorithm to a many-electron system yields a totally symmetric 
many-boson ground state of no physical interest. The fixed-node approximation introduces a 
trial many-electron ground-state wavefunction, ^^(R), and forbids walker moves that cause 
Tr to change sign. As long as Tt is properly antisymmetric, this is sufficient to ensure that a 
fermionic solution is obtained. It may be showrP^that energies calculated within the fixed- 
node approximation are variational: the result is greater than or equal to the many-fermion 
ground-state energy and tends to the exact energy as the (3A"e —l)-dimensional nodal surface 
on which = 0 tends to the ground-state nodal surface. The fixed-node approximation is 
required for DMC simulations of large systems but is the only fundamental limitation of the 
method. Other approximations, such as the use of a finite time-step or the representation of 
ions by pseudopotentials, can be made negligible or avoided given sufficient computer time. 
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Fixed-node DMC energies are in most cases comparable in accuracy to energies calculated 
using the CCSD(T) method,^ which is often known as the “gold standard” of quantum 
chemistry. 

The diffusion/branching simulation described above is unstable in practice because the 
potential energy V (R) diverges whenever electrons approach nuclei or each other, leading 
to uncontrollable branching. This problem can be overcome using an importance-sampling 
technique. The imaginary-time Schrodinger equation is re-expressed in terms of the quantity 
/(R, r) = Tr(R)T(R, r) to obtain 

r) - Vr . |v(R)/(R, t)] - |i54R) - Bt|/(R, r), (3) 

where Vr = (Vn, Vra, ■ ■ ■, Vr^^) is the 3iVe-dimensional gradient operator, Vi, = Vr ■ 
Vr is the corresponding Laplacian, v(R) = VrIu |\['r(R)| is the 3iVe-dimensional drift 
velocity vector, and i?L(R) = (l/T' 7 ’(R))RTr(R) is the local energy. The importance- 
sampled imaginary-time Schrodinger equation describes a diffusion process similar to that 
discussed above, except that the walkers now drift with velocity v(R) as well as diffusing and 
branching. The branching rate is determined by the shifted local energy El(R.) — Et instead 
of the shifted potential energy ld(R) — Erp. If the trial function is a good approximation to 
the ground state, the local energy is a smooth function of R and the numerical difficulties 
caused by divergences in ld(R) are avoided. The fixed-node approximation is imposed by 
rejecting walker moves that change the sign of ^^(R). 

Our DMC simulations used the CASINO QMC cod^^S ^ trial function of Slater- 
Jastrow (SJ) form, 

Tr(R) = exp[J(R)] det[?/>„(rj)] det[^/’„(rj)], (4) 

T I 

where rj is the position of the i’th spin-up electron, is the position of the j’th spin-down 
electron, exp[J(R)] is the Jastrow factor, and det['0n(i'I)] and det['0„(rj)] are Slater deter¬ 
minants of spin-up and spin-down one-electron orbitals. These orbitals were obtained from 
DFT calculations using the plane-wave-based Quantum Espresso code.^^ A norm-conserving 
pseudopotential constructed within DFT using the Perdew-Zunger parameterization of the 
local density approximatiorP^ was employed. We chose a very large basis-set cut-off of 
300 Ry to guarantee converge to the complete basis-set limit.^ The one-electron orbitals, 
originally expressed as linear combinations of plane waves, were transformed into a blip poly¬ 
nomial basis for efficiency.^ The Jastrow function J(R) consisted of polynomial one-body 
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electron-nucleus (en) and two-body electron-electron (ee) terms, the parameters of which 
were optimized by variance minimization at the variational Monte Carlo (VMC) level.^^^^ 

B. Finite-size errors and correction methods 

QMC FS errors are conventionally separated into one-body and many-body contributions. 
One-body (independent-particle) errors arise from the non-interacting kinetic, potential and 
Hartree energies and include shell-hlling effects. These errors can be removed by twist 
averaging.!!^ Many-body errors arise from the effects of exchange and correlation on the 
Coulomb and kinetic energies and are not removed by twist averaging. As explained in the 
introduction, various techniques may be used to reduce or cancel these errors, but none 
is entirely successful and care is required. The oldest approach is extrapolation, which 
remains useful. The use of the modihed periodic Coulomb interactiorP^^^^^ reduces the 
Coulomb errors but not the kinetic energy errors and must therefore be combined with other 
techniques. The LDA-based Kwee-Zhang-Krakauer (KZK) approach!^ applies corrections 
obtained from DFT calculations carried out using a modihed exchange-correlation functional 
explicitly designed to mimic the DMC many-body errors. 

To remove single-particle errors and eliminate shell effects in the kinetic energies of metal¬ 
lic systems, we use twist-averaged boundary conditions.^ A twist is imposed by insisting 
that the many-electron wavefunction obeys the Bloch boundary condition 

Tk,(ri,..., Fj + L,..., rjvJ = exp(ik, ■ L)Tk,(ri, ..., r^, ..., rjvJ (5) 

for all electrons i, where L is any simulation-cell lattice vector. Expectation values of 
observables are obtained by averaging over twist vectors k^ uniformly distributed over the 
simulation-cell Brillouin zone: 

= ( 6 ) 

-'''twist , 

ks 

The twists can be chosen from a uniform Monkhorst-Pack grid,^ preferably offset from F, 
or can be chosen randomly, as in this work. The number of twists should be as large as 
computational resources allow. 

As the twist k^ varies, the energies of some of the one-electron states appearing in the 
Slater determinants may cross the Fermi level. In the canonical approach to twist averaging. 
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the electron number is kept constant and the Fermi level is allowed to vary with twist. This 
makes the twist-averaged total energy slightly too large,the bias reduces as the 
simulation-cell increases in size and is normally negligible. In the grand canonical approach 
to twist averaging, the Fermi energy is hxed and the electron number is allowed to vary with 
ks. As was demonstrated in Ref. [101 energies obtained using grand-canonical twist averaging 
exhibit much larger fluctuations than energies obtained using canonical twist averaging with 
the same number of twists. Furthermore, these fluctuations die away very slowly as the 
number of twists is increased. To save effort, applications of QMC to real systems normally 
use the smallest number of twists possible, so canonical twist averaging is preferable despite 
the fact that it suffers from small systematic errors. This work only considers twist averaging 
within the canonical ensemble. 

In metallic systems, even when substantial computational resources are expended, the FS 
errors due to incomplete twist averaging are substantial. We therefore dehne an incomplete- 
twist-averaging correction as follows: 

Abz = E^^^{oo)-E^^^{L), ( 7 ) 

where E^^^{oo) is the DFT energy computed using a fully converged fc-point mesh and 

E?av(.L) = ^ E k.) (8) 

ks 

is the twist-averaged DFT energy obtained using the same simulation cell and set of twists as 
the DMC simulation. The incomplete-twist-averaging correction tends to zero as the DMC 
twist averaging is improved and works well if the independent-particle hnite-size errors are 
well approximated by their DFT equivalents. In practice, this approach allows accurate 
results to be obtained with surprisingly small sets of DMC twists, even in metals. 

We analyze three different methods for correcting the many-body FS errors in DMC 
results. Two of these use the structure-factor-based corrections proposed by Chiesa et 
The hrst employs the standard Ewald form of the periodic Coulomb interaction and Chiesa 
corrections for both the kinetic and potential energies; the second uses the MPCP^*^^*^ to 
deal with the Coulomb errors and a Chiesa correction for the kinetic energy only. Results 
obtained with both of these methods are expected to be similar in quality.^ 

The third FS-correction method considered here is the KZK approach,!^ which uses a 
system-size-dependent local density approximation htted to the results of DMC simulations 









of finite cubic simulation cells of uniform electron gas. DFT energies calculated using the 
KZK functional incorporate DMC FS errors within an approximation analogous to the LDA. 
To estimate the FS error in the DMC total energy of a given simulation cell, the DFT total 
energy of exactly the same simulation cell is calculated using the KZK functional. The 
difference between this value and the DFT energy of an inhnite simulation cell calculated 
using the standard LDA provides an estimate of the DMC FS error. The KZK functional 
was not originally combined with twist averaging, but the combination is easy to implement 
(see below) and very successful. 

In the following we explain how to combine twist averaging and KZK corrections. In 
general, any FS correction can be written as 

AF;^"(L) = Af|(L) + A^^^(L), (9) 

where Affi(L) includes contributions from the Hartree energy, the electron-nuclear Coulomb 
interaction energy, and the one-body component of the kinetic energy, while A^^(L) is a 
many-body term that includes contributions from the exchange-correlation energy and the 
many-body part of the kinetic energy. More precisely, Affi (L) may be dehned as that part 
of the total hnite-size error that is also present in a DFT calculation for the same simulation 
cell and can be corrected using DFT results. 

In their original paper,^ Kwee, Zhang and Krakauer considered the hnite-size errors 
affecting a QMC simulation carried out in a supercell of L x L x L primitive unit cells with 
twist ks = 0. The corresponding one-particle hnite-size error is 

A(i{L) = E^^^{oo)-E^^^{L), ( 10 ) 

where E^^'^{oo) is the DFT energy obtained using a fully converged fc-point mesh and 
E^^'^{L) is the F-point DFT energy of the supercell. This, of course, can be calculated 
using an L X L X L Monkhorst-Pack grid of k points in the primitive Brillouin zone. Since 
ks = 0, the Monkhorst-Pack grid includes the origin. The KZK approximation to the 
many-body hnite-size error is 

A^J^iL)^E^^^{L)-E^^^iLfi ( 11 ) 

where E^^^{L) is the P-point DFT energy of the supercell computed using the KZK func¬ 
tional instead of the standard LDA. The KZK approximation to the total hnite-size error 
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is 

AE^^(L) = Afl(L) + A^^^(L) ^ - E^^^{L). (12) 

If twist averaging is used, the DMC energy becomes a function of the twist k*, which lies in 
the small Brillouin zone corresponding to the simulation supercell. The Slater determinants 
appearing in the twisted trial wavefunction are built using orbitals from a Monkhorst-Pack 
grid oi L X L X L k points within the larger primitive Brillouin zone, offset by from 
the origin. In a DFT context, carrying out a supercell calculation at a non-zero twist k^ is 
equivalent to approximating the integration over the primitive Brillouin zone by a quadrature 
over this offset grid of k points. 

To help analyze the FS errors, we write the DMC ground-state energy of the inhnite 
simulation cell as 

^ E^^y^iL) + AE^UL): (13) 

where 

(14) 

twist I 

ks 

is the twist-averaged DMC energy of the L x L x L simulation cell and AE^^y{L) is the 
required FS correction. In the spirit of KZK, this is approximated using the formula 

AE^^L) ^ E^^^ioo) - E^^^iL), (15) 

where E^^'^{oo) is the DFT energy computed within the LDA using a fully converged 
/c-point mesh (which in this work means 28 x 28 x 28) and 

e¥IP'{L) k.) (16) 

-'''twist , 

ks 

is the twist-averaged KZK energy for the supercell, computed using the same set of Atwist 
twists employed in the DMC simulations. The FS correction, AE^^y{L), accounts both for 
the many-body FS errors and for any one-body FS errors not removed by the limited twist 
averaging employed in the DMC simulations. 

The use of twisted boundary conditions requires the use of complex trial wavefunctions 
and increases the computational cost a little because complex arithmetic is slower than real 
arithmetic. In the VMC and wavefunction optimization algorithms, since the expectation 
values of Hermitian operators must be real, only the real parts of the local-energy components 
need to be calculated and collected. The run-time and programming-time costs of twist 
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averaging are therefore small. The use of a complex trial function in DMC requires the 
replacement of the hxed-node approximation, in which the DMC wavefunction is constrained 
to have the same sign as the trial wavefunction, by the hxed-phase approximation,l2Slin which 
the DMC and trial wavefunctions are constrained to have the same phase. In practice, 
however, the hxed-node and hxed-phase algorithms are very similar and little extra coding 
is required: the real part of the drift vector is used when proposing trial electron moves; it 
is neither necessary nor possible to reject node-crossing electron moves; and, as in VMC, 
only the real parts of the local energies are gathered. Another important issue in twist 
averaging is the Jastrow factor. This work uses the same optimized Jastrow for each twist 
vector, as we found that re-optimizing the Jastrow factor at every twist provides negligible 
improvements in the hnal results.^ We note, hnally, that the VMC or DMC runs at each 
twist can be relatively short and need not be fully converged. The idea is that we collect 
enough data to achieve an acceptable error bar when the data are averaged over all twist 
vectors. If a normal run without twist averaging takes N moves to arrive at an acceptable 
error bar, each twist angle need only be run for around V/Wwist moves. 


III. RESULTS AND DISCUSSION 


This section presents DMC results obtained using the three different FS-correction meth¬ 


ods explained in the Sec. IIB As simple example metals, we have studied lithium (Li) and 


aluminum (Al), with one and three valence electrons, respectively. The frozen ionic cores 
are represented by non-local norm-conserving LDA pseudopotentials. The KZK functional 
is essentially an LDA, so the use of LDA pseudopotentials allows us to obtain a consistent 
comparison of all three hnite-size-correction approaches considered. There is evidence that 
Hartree-Fock pseudopotentials may produce more accurate results than LDA pseudopoten¬ 
tials when used in DMC simulations, but since our aim is to investigate FS errors, and since 
these are almost independent of pseudopotential, no advantage would be gained by using 
another pseudopotential type. To check the accuracy and convergence of the DMC energy as 
a function of the number of atoms N in the simulation cell, we have performed calculations 
for a range of different values of N (and thus also different values of L). 

We hrst investigate the effects of applying Chiesa and MFC corrections to the results of F- 
point DMC simulations with no twist averaging. Table |T]shows the F-point DMC results for 
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TABLE I. DMC energies of metallic lithium for different numbers N of atoms in the simulation 
cell. Every energy appearing in the table is the outcome of a single E-point DMC simulation for 
the appropriate simulation cell. Results obtained by applying Chiesa kinetic (AKE) and potential 
(APE) hnite-size corrections to DMC energies calculated using the Ewald interaction agree well 
with results obtained by applying only the Chiesa kinetic correction to DMC energies calculated 
using the modihed periodic Coulomb interaction. Despite the application of many-body hnite-size 
corrections, the calculated energy depends strongly on the size of the simulation cell. This indicates 
that the single-particle hnite-size errors are large. Energies are in eV per atom. 


N Ewald Ewald -|- AKE Ewald -1- APE Ewald -1- AKE -I- APE 

MPC MPC -b AKE 

32 -6.94307(3) -6.92246(3) 

-6.87450(3) 

-6.85389(3) 

-6.86219(4) -6.84157(4) 

48 -6.97920(7) -6.97852(7) 

-6.91729(7) 

-6.91661(7) 

-6.91661(6) -6.91593(6) 

72 -6.98253(3) -6.97981(3) 

-6.94471(3) 

-6.94205(3) 

-6.94062(3) -6.93790(3) 

96 -6.97063(5) -6.96845(5) 

-6.94015(5) 

-6.93804(5) 

-6.94382(5) -6.94171(5) 

144 -6.89491(2) -6.89089(2) 

-6.87634(2) 

-6.87232(2) 

-6.87640(2) -6.87246(2) 


Li obtained using various different FS correction methods. As expected,^ results (denoted 
Ewald -|- AKE -|- APE) obtained by adding Chiesa kinetic and potential energy corrections 
to the DMC energy calculated using the Ewald interaction are in good agreement with results 
(denoted MPC -|- AKE) obtained using the modihed periodic Coulomb interaction with a 
Chiesa correction for the kinetic energy only. Note that the higher-order kinetic energy 
corrections dehned according to equation (55) in Ref. dUl are included in AKE. Because of 
the lack of twist averaging, the single-body FS errors are large and the calculated ground- 
state energies depend strongly on the size of the simulation cell. The choice of the F-point, 
ks = 0, maintains the symmetry of the system but usually increases shell-hlling effects in 
metallic systems, making the independent-particle FS errors even worse. It is clear, however, 
that no calculation carried out at a single twist vector will yield satisfactorily accurate results. 
The use of twist averaging is essential in metals. 

Tables |TT] and present twist-averaged DMC results for Li and Al, respectively, again 
corrected using the Chiesa and MPC approaches. The integration over the simulation- 
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cell Brillouin zone that produces the twist-averaged energy would completely remove the 
single-particle FS errors if carried out exactly, but in practice the integration has to be 
approximated as a summation over a hnite set of twists. The summation tends to the 
integral as either the number N of atoms in the simulation cell or the number iVtwist of twists 
used tends to inhnity, but is far from perfect in practice. Checking the convergence with 
respect to both N and A^twist is important, since fully twist-averaged calculations for hnite 
simulation cells still contain many-body FS errors and do not necessarily provide accurate 
results. The main reason is that the fully twist-averaged exchange-correlation energy still 
depends on N, even though the one-electron part of the fully twist-averaged kinetic energy 
does not. Comparing the hnite-size-corrected and twist-averaged DMC energies from Tables 
[IT] with the F-point energies from Table [I] shows that twist averaging much reduces the hnite- 
size errors and allows accurate results to be obtained using much smaller simulation cells. 
In this particular case (but not in general) it also produces a lower ground state energy. 

Our twist-averaged DMC results were obtained using 24 randomly sampled twists (values 
of ks) in the simulation-cell Brillouin zone. Two other practical sampling schemes exist. 
One is to use a uniform Monkhorst-Pack gricP^ of points centered on the F-point of 
the simulation-cell Brillouin zone, and the other is to use a uniform grid centered on the 
Baldereschi poinlP^ of the simulation-cell Brillouin zone. As either the number of twists or 
the size of the simulation cell tends to inhnity, all three twist-averaging methods yield the 
same results. 

In applications of DMC to real systems using computers routinely available to researchers, 
it is rarely possible to treat very large simulation cells or numbers of twists. Restricting the 
number of twists is particularly problematic in metallic systems, where the Fermi surface 
discontinuity makes the integrand (for example, the total kinetic energy) a discontinuous 
function of k^ at zero temperature. The convergence with system size and number of twists 
is therefore much slower for metals than for insulators. Hartree-Fock calculations for the 
uniform electron ga^ show that energies obtained using Baldereschi twist averaging con¬ 
verge faster than energies obtained using random twist averaging at very large numbers of 
twists (although it could be argued that the choice of the Baldereschi point introduces a 
systematic bias into the unconverged results), but that both methods converge slowly. Here 
we show that the use of the incomplete-twist-averaging correction dehned in Eq. ([^ allows 
well-converged results to be obtained with very small numbers of twists. The choice between 
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TABLE II. DMC energies of lithium for different numbers N of atoms in the simulation cell. 
Every energy appearing in the table is an average of the outcomes of 24 separate DMC simulations 
with different randomly chosen twists. Results obtained by applying Chiesa kinetic (AKE) and 
potential (APE) finite-size corrections to DMC energies calculated using the Ewald interaction 
agree well with results obtained by applying only the Chiesa kinetic correction to DMC energies 
calculated using the modified periodic Coulomb interaction. The use of twist averaging has much 
reduced the independent-particle finite-size errors observed in TableEnergies are in eV per atom. 


N Ewald 

Ewald -|- AKE Ewald -|- APE Ewald -|- AKE -|- APE 

MPC 

MPC -7 AKE 

32 -7.0008(4) 

-6.9802(4) 

-6.9334(4) 

-6.9128(4) 

-6.9301(6) 

-6.9095(6) 

48 -6.9708(4) 

-6.9699(4) 

-6.9100(4) 

-6.9092(4) 

-6.9089(5) 

-6.9080(5) 

72 -6.9590(3) 

-6.9563(3) 

-6.9207(3) 

-6.9180(3) 

-6.9203(2) 

-6.9177(2) 

96 -6.9489(2) 

-6.9467(2) 

-6.9184(2) 

-6.9162(2) 

-6.9207(2) 

-6.9176(2) 

144 -6.9291(2) 

-6.9251(2) 

-6.9113(2) 

-6.9073(2) 

-6.9118(2) 

-6.9078(2) 


Baldereschi or random twist sampling is then unimportant. 

Tables [IV] and |V] compare twist-averaged DMC results obtained using the ABZ incomplete- 
twist-averaging FS correction and several different many-body FS-correction methods. The 
convergence of the FS-corrected DMC energies with system size N is excellent and there 
is no difficulty in reaching an accuracy of a few meV per atom. As before, our DMC 
energies are averages over 24 randomly chosen twists in the simulation-cell Brillouin zone, 
corresponding to 24 randomly translated L x L x L Monkhorst-Pack fc-point meshes in the 
primitive Brillouin zone. We also carried out twist-averaged DMC calculations using 36 
twists; the change in the total energy was less than 1.5 meV/atom for Li and less than 2.7 
meV/atom for Ah This shows that twist-averaged DMC energies including ABZ corrections 
converge very rapidly as the number of twists is increased. 

Figure shows how the FS-corrected DMC energies of metallic Li depend on system 
size, allowing an easy comparison of the three different many-body FS-correction methods 
considered in this work. All DMC energies are averaged over 24 randomly chosen twists 
and include ABZ corrections. Red squares indicate DMC results calculated using the Ewald 
interaction with Chiesa corrections for the kinetic and potential energies (Ewald -|- AKE 
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TABLE III. DMC energies of aluminum for different numbers N of atoms in the simulation cell. 
Every energy appearing in the table is an average of the outcomes of 24 separate DMC simulations 
with different randomly chosen twists. Results obtained by applying Chiesa kinetic (AKE) and 
potential (APE) finite-size corrections to DMC energies calculated using the Ewald interaction 
agree well with results obtained by applying only the Chiesa kinetic correction to DMC energies 
calculated using the modified periodic Coulomb interaction. Energies are in eV per atom. 


N Ewald Ewald -|- AKE Ewald -|- APE Ewald -|- AKE -|- APE 

MPC 

MPC -h AKE 

24 -56.175(2) 

-56.170(2) 

-55.932(2) 

-55.927(2) 

-56.095(2) 

-56.091(2) 

32 -56.203(1) 

-56.152(1) 

-56.074(1) 

-56.022(1) 

-56.0585(2) 

-56.003(2) 

48 -56.155(1) 

-56.138(1) 

-56.048(1) 

-56.031(1) 

-56.058(2) 

-56.041(2) 

72 -56.0922(7) 

-56.0808(7) 

-56.0212(7) 

-56.0098(7) 

-56.0226(8) 

-56.0112(8) 


TABLE IV. DMC energies of lithium for different numbers N of atoms in the simulation cell. 
Every DMC energy is an average of the outcomes of 24 separate DMC simulations with different 
randomly chosen twists. The incomplete-twist-averaging finite-size correction ABZ is included in 
all energies. Results obtained by applying the Chiesa kinetic (AKE) and potential (APE) finite- 
size corrections to DMC energies calculated using the Ewald interaction agree well with results 
obtained by applying the Chiesa kinetic correction to DMC energies calculated using the modified 
peiodic Coulomb interaction. Results obtained using the twist-averaged KZK method, which also 
include an equivalent of the ABZ correction, are also in good agreement. Energies are in eV per 
atom. 


N Ewald AKE APE -h ABZ MPC AKE ABZ TAV-KZK 


32 

-6.9299(4) 

-6.9267(4) 

-6.9126(4) 

48 

-6.9183(4) 

-6.9172(4) 

-6.9095(4) 

72 

-6.9186(3) 

-6.9182(3) 

-6.9126(3) 

96 

-6.9149(2) 

-6.9163(2) 

-6.9123(2) 

144 

-6.9142(2) 

-6.9147(2) 

-6.9125(2) 
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TABLE V. DMC energies of aluminum for different numbers N of atoms in the simulation cell. 
Every DMC energy is an average of the outcomes of 24 separate DMC simulations with different 
randomly chosen twists. The incomplete-twist-averaging finite-size correction ABZ is included in 
all results. Results obtained by applying the Chiesa kinetic (AKE) and potential (APE) finite-size 
corrections to DMC energies calculated using the Ewald interaction agree well with results obtained 
by applying the Chiesa kinetic correction to DMC energies calculated using the modified peiodic 
Coulomb interaction. Results obtained using the twist-averaged KZK method, which also include 
an equivalent of the ABZ correction, are also in good agreement. Energies are in eV per atom. 


N 

Ewald -F AKE -h APE -h ABZ MPC AKE -h ABZ 

TAV-KZK 

24 

-56.019(2) 

-56.183(2) 

-56.025(2) 

32 

-56.088(1) 

-56.069(2) 

-56.088(1) 

48 

-56.078(1) 

-56.088(2) 

-56.081(1) 

72 

-56.0607(7) 

-56.0621(8) 

-56.0626(7) 


-|- APE -|- ABZ); green circles indicate DMC resnlts obtained nsing the modified periodic 
Conlomb interaction with Chiesa corrections for the kinetic energy only (MPC -|- AKE 
-|- ABZ); and blue circles indicate KZK-corrected DMC results, again incorporating ABZ 
corrections. Even for the smallest simulation cell considered, with just 32 atoms, the errors 
in the Ewald -|- AKE -|- APE -|- ABZ and MPC -|- AKE -|- ABZ total energies are only 
15.7 and 12 meV/atom, respectively. The errors in Ewald DMC energies corrected using 
the KZK scheme are even smaller, at approximately 3 meV/atom. 

We emphasize the importance of the success of the twist-averaged KZK method from 
a practical point of view. It is known, for example, that the cheap and widely used DFT 
approach often fails to provide accurate enough resultd^ to understand the behavior of 
materials at high pressure. Therefore, to study the very interesting phase diagram of Li,^ 
it will be necessary to perform full many-body computations, most likely using DMC. The 
drawback is that DMC calculations are typically at least a thousand times more expensive 
than DFT calculations. The twist-averaged KZK approach allows one to investigate a large 
number of possible crystal structures and construct the Li phase diagram whilst keeping the 
cost of the DMC simulations within reasonable bounds. 
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FIG. 1. Total DMC energies of Li as function of the number of particles in the simulation 
cell. All energies are averaged over 24 randomly chosen twists and include ABZ corrections. Red 
squares indicate DMC results calculated using the Ewald interaction with Chiesa corrections for 
the kinetic and potential energies (Ewald + AKE + APE + ABZ); green circles indicate DMC 
results obtained using the modified periodic Coulomb interaction with Chiesa corrections for the 
kinetic energy only (MPC + AKE + ABZ); and blue circles indicate KZK-corrected DMC results, 
again incorporating ABZ corrections. 

Figure shows how the FS-corrected DMC energy of metallic A1 depends on the number 
of atoms in the simulation cell. All DMC energies are averaged over 24 randomly chosen 
twists and include ABZ corrections. Red squares indicate FS-corrected results obtained 
using the Ewald interaction with Chiesa corrections for the kinetic and potential energies 
(Ewald -|- AKE -|- APE -|- ABZ); green circles indicate DMC results obtained using the 
modihed periodic Coulomb interaction (MPC) with Chiesa corrections for the kinetic energy 
only (MPC -|- AKE -|- ABZ); and blue circles indicate KZK-corrected DMC results, again 
incorporating ABZ corrections. For all simulation-cell sizes, the Ewald -|- AKE -|- APE -|- 
ABZ and KZK results are in almost perfect agreement. The difference between the MPC 
-|- AKE -|- ABZ energies and those obtained using the other two FS-correction methods is 
about 7 meV/atom for a cell containing just 48 atoms and decreases rapidly with increasing 
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simulation-cell size. 

The results of this section have shown that the addition of an incomplete twist-averaging 
correction, ABZ, allows accurate results to be obtained with a remarkably small number of 
twists, even for metals. In most of the cases studied (except for the 96-atom Li simulation 
cell), the ABZ correction lowers the total DMC energy. The values of ABZ for Li simulation 
cells containing 32, 48, 72, 96, and 144 atoms are -0.0171, -0.0091, -0.0006, -1-0.0013, and 
-0.0069 eV, respectively. In the case of A1 simulation cells containing 24, 32, 48, and 72 
atoms, the values of ABZ are -0.092, -0.066, -0.047, and -0.050 eV, respectively. A recent 
paper by Shulenberger and MattssorP^ provided accurate benchmark DMC results for a 
wide range of different bulk materials. They required 216 and 64 twists to obtain converged 
results for Li and A1 supercells containing 28 and 108 atoms, respectively. Because we use 
incomplete twist-averaging corrections, ABZ, we are able to obtain similarly accurate results 
with considerably fewer twists. 


IV. CONCLUSION 

We have systematically analyzed and compared the various schemes that have been pro¬ 
posed for correcting FS errors in QMC simulations of real metallic systems. We have ex¬ 
plained how to combine the use of twist-averaged boundary conditions with the KZK func¬ 
tional and shown the value of incomplete-twist-averaging corrections based on DFT. The 
reassuring news is that all of the commonly used approaches work well. 

We believe that the use of DFT-based incomplete-twist-averaging corrections will have 
an important role to play in DMC simulations of real metallic systems. The reliance on 
DFT could be considered a drawback, but it is important to bear in mind that any valid FS- 
correction method must yield the same total energy in the limit of a large enough simulation 
cell. The important question is not whether the unattainable limiting value is correct but 
how rapidly it is approached. The use of incomplete-twist-averaging corrections significantly 
improves this convergence. Furthermore, energies calculated using the twist-averaged KZK 
scheme (which implicitly incorporates a ABZ incomplete-twist-averaging correction) often 
settle down to a system-size-independent constant more quickly than energies calculated 
using other methods incorporating ABZ corrections. 

We believe that this paper will provide a useful guide and benchmark for researchers 
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1/N 

FIG. 2. Total DMC energies of A1 as function of the number of particles in simulation cell. All 
energies are averaged over 24 randomly chosen twists and include ABZ corrections. Red sqnares 
indicate DMC results calculated using the Ewald interaction with Chiesa corrections for the kinetic 
and potential energies (Ewald + AKE + APE + ABZ); green circles indicate DMC results obtained 
using the modified periodic Coulomb interaction with Chiesa corrections for the kinetic energy only 
(MPC + AKE + ABZ); and bine circles indicate KZK-corrected DMC results, again incorporating 
ABZ corrections. 

using QMC and other many-body electronic structure methods such as CCSD(T) and the 
GW approximation to study metallic systems. 
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